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SUMMARY 


The ION computer code is designed to calculate charge exchange 
ion densities, electric potentials, plasma temperatures, and current 
densities external to a neutralized ion engine in R-Z geometry. The 
present version assumes the beam ion current and density to be known 
and specified, and the neutralizing electrons to originate from a 
hot-wire ring surrounding the beam orifice. The plasma is treated as 
being resistive, with an electron relaxation time comparable to the 
plasma frequency. Together with the thermal and electrical boundary 
conditions described below and other straightforward engine 
parameters, these assumptions suffice to determine the required 
quantities. 

The ION code, written in ASCII FORTRAN for UNI VAC 1100 series 
computers, is designed to be run interactively, although it can also 
be run in batch mode. The input is free-format, and the output is 
mainly graphical, using the machine-independent graphics developed for 
the NASCAP code. The executive routine calls the code's major 
subroutines in user-specified order, and the code allows great 
latitude for restart and parameter change. 
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1. INTRODUCTION 


The ION computer code is designed to calculate charge exchange 
ion densities, electric potentials, plasma temperatures, and current 
densities external to a neutralized ion engine in R-Z geometry. The 
present version assumes the beam ion current and density to be known 
and specified, and the neutralizing electrons to originate from a 
hot-wire ring surrounding the beam orifice. The plasma is treated as 
being resistive, with an electron relaxation time comparable to the 
plasma frequency. Together with the thermal and electrical boundary 
conditions described below and other straightforward engine 
parameters, these assumptions suffice to determine the required 
quantities. 

The ION code, written in ASCII FORTRAN for UNIVAC 1100 series 
computers, is designed to be run interactively, although it can also 
be run in batch mode. The input is free-format, and the output is 
mainly graphical, using the machine-independent graphics developed for 
the NASCAP code. The executive routine calls the code's major 
subroutines in user-specified order, and the code allows great 
latitude for restart and parameter change. 

In Chapter II of this report we outline the theoretical 
treatment employed by the ION code. Chapter III discusses the use of 
the code, describing the input parameters and the functions of the 
major subroutines. Chapter IV presents some sample results, and 
Chapter V contains a sample run. 
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2. THEORY 


The major theoretical components to the ION computer code have 
been previously presented in the open literature. This Chapter will 
outline and review the code theory. For convenience, the published 
detailed descriptions are reproduced in Appendices A and B. 

A 30 cm diameter mercury thruster with 2 A of beam current and 

25 mA of charge exchange ion current can be expected to produce charge 

14 ^3 

exchange ion densities at the beam edge on the order of 10 m 
(Appendix B and references therein). Since plasma densities in low 
earth orbit are on the order of 10^^ - 10 m , the charge 
exchange plasma will dominate the local environment. Such a plasma 
will be characterized by long collision lengths and ( ~lo\) and 

_3 

short Debye lengths ( -10 m). 

While this plasma is not collision dominated, the short Debye 
length implies that small scale fluctuating fields may create an 
effective collision frequency, that will be effective in 

randomizing electron velocities on the scale of typical problem 
dimensions ( ~lm). This suggests that except within a few Debye 
lengths of particle sinks and sources (or the sheaths of high voltage 
surfaces) we can expect the plasma to be quasineutral. Quasineutrality 
is a basic assumption throughout ION. Appendix A describes the 
theoretical and experimental evidence for a barometric law correlation 
between potential d, density p, and temperature e: 

d(r)=ejin(^h . 

^0 

Also described in Appendix A are experiments where a barometric law 
did not adequately explain the observations, suggesting the need for a 
comprehensive electron transport model. ION utilizes a barometric law 
assumption only as an initial approximation. In Appendix A it is 


4 


demonstrated that measured neutralizer electron temperatures and 
potential variations can be qualitatively explained in terms of the 
anomalous resistivity of the thruster generated plasma to the flow of 
electrons if we invoke an effective collision frequency comparable to 
the electron plasma frequency. 

The ION code is constructed from four main theoretical models: 

An ion beam model, a free expansion model for neutrals, a hydrodynamic 
charge exchange ion model, and a fluid model for the neutralizer 
electrons. Of these four, the first two are analytic and require very 
little computing while the second two are numerical. The neutral 
expansion model is described in Appendix C, the hydrodynamic charge 
exchange ion model in Appendix B, and the electron fluid model in 
Appendix A. The code operates in R-Z geometry on a discrete rectangu- 
lar mesh of nodal points. At present, the ion beam is modeled as 
columnar with no spreading. Depletion due to charge exchange is 
ignored. The beam originates at z = 0, is centered at r = 0, and has 
either a quadratic or Gaussian profile. (See Chapter III). The 
neutral efflux is described analytically by a free expansion model 
assuming a uniform emission over the thruster throat with a Maxwellian 
velocity spread. The charge exchange ions are generated within the 
beam with zero initial velocities. The generation rate is equal to 
the product of the local beam flux, neutral density, and the charge 
exchange cross-section. 

We use a hydrodynamic description for both the charge exchange 
ions and the neutralizer electrons. However, for each case there are 
differing approximations that can be made in solving the hydrodynamic 
equations. Therefore, purely for the purpose of discussion, the charge 
exchange model is referred to as hydrodynamic, and the electron model 
is labeled fluid. 

The hydrodynamic model of the charge-exchange ions, described in 
Appendix B, is based upon a time stepped solution of the mass and 
momentum equations in finite difference form. A barometric law 
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relationship is assumed between density, potential and temperature. 

The mechanics of the finite difference formulations are given in 
detail in Appendix B. The major simplifying assumption in this model 
is that charge exchange ions are quite cold so that pressure terms may 
be neglected. This means that the inertial term in the momentum 
equation ( v *p v v) in Eq. (3) of Appendix B) forces the equations to 
be parabolic in nature. It is this parabolic nature that necessi- 
tates a time-stepped finite difference approach. The initial and 
boundary conditions for this model are discussed in Chapter III. 

The electron fluid calculation follows after the charge exchange 
ion calculation using the ion densities as input. For the electron 
model, we make the assumption that the plasma offers a significant 
resistance to the flow of electrons (due to the effective collision 
frequency). It then follows that the electron drift velocities will be 
small compared to thermal velocities, which allows the inertial term 
in the momentum and energy equations to be dropped. Thus, the 
conservation equations for mass, momentum and energy may be solved, 
along with the ideal gas equation of state, in closed elliptic form. 

A finite element approach is used to solve directly for the steady 
state. This model is discussed in Appendix A. 

The switching between finite difference and finite element 
approaches may seem cumbersome, but the conversions are quite efficient 
compared to the difficulties encountered otherwise. Finally, it 
should be mentioned that since there is currently no iteration between 
the ion and electron calculations, the resulting parameters will not 
be entirely self-consistent. This limitation does not appear to be 
severe in the test cases modeled so far, but a truly predictive model 
should be free from such limitations. This will be a subject for 
future work. 
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3. CODE USAGE 


A. OVERVIEW 

The ION code is invoked by the statement 
0 XQT absol uteel ementname 

It assigns files 9 and 21 as restart files, which may be aliased 
(0USE) to user permanent files, or later copied to permanent files. 

It also assigns and writes printed output on files (by default) 19 and 
20, and uses files 11 and 12 as scratch files. Graphical output is 
generated by writing NASCAP-type pseudo plot calls on file 2 and 
interfacing to the user's graphics library with the standard NASCAP 
postprocessor (PLOTREAD). 

After invoking the ION code, the user's first task is to enter 
the problem parameters (Table 3.1). Acceptable sets of input are 
shown in Figures 3.1. After encountering an 'END' card in the 
parameter input, a summary is printed (Figure 3.2). 

Subsequent major subroutine calls are made at user request. 

These calls are listed in Table 3.2. For a new problem the first call 
will be to subroutine CHEX, which determines the steady state 
charge-exchange ion density. (This calculation usually converges in 
the default of 300 steps, but may take several minutes of computer 
time. For batch runs, a larger value of 500 steps is recommended.) 
Next subroutine POTENT is called to perform the calculation of 
temperatures and potentials. Entry points RECHEX and RESUME are 
available to restart a non-converged calculation of CHEX or POTENT, 
respectively, and entry point NEWPOT may be used instead of POTENT to 
perform a calculation with altered thermal boundary conditions or 
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TABLE 3.1. PROBLEM PARAMETERS 


Grid Size 


KEYWORD 

DEFAULT 

MEANING 

NR 

25 (max) 

Number of radial coordi- 
nate values. 

NZ 

25 (max) 

Number of Z-coordinate 
values. 

RHAX 

1. 

Radius of computational 
space, (meters) 

ZMAX 

1. 

Length of computational 


Ion Beam Characteristics 


space, (meters) 


CURRent 

I. 

Beam current, (amperes) 

RADIUS 

I. 

Beam radius on entry, (m) 

ENERgy 

1000. 

Beam Ion energy (eV). 

MASS 

200.59 

Ion mass In amu. 

PROFIIe 

GAUS 

Beam profile acceptable 
values are 'GAUS' and 



'QUAD' . 

SPREAdangle 

0. 

Half-angle character- 


Izing beam spread. (Non- 
zero spread angle not 
Implemented.) 


Charge Exchange Ion Production (Physical) 
SIGMA 6x10-1 


6x10-19 Charge exchange produc- 
tion cross-section (m‘) 


FLOW 2. 

(variable CNEUT) 

TEMPerature 1. 


Fuel flow In ampere 
equivalent 

Neutral temperature (eV). 


keyword default 

Charge Exchange Ion Production (Computational) 
ECHA 10. 


MEANING 


VDTR 0.2 

(variable VDTBYR) 

eta .05 

RHOMIn l.El 

NSTEPs 300 

Thermal and Current Boundary Conditions 
THETal 1. 

THERbc 'I SOT 


RNEUtr 

NlTEr 


Output Destinations 


Characteristic accelera- 
ting voltage for charge 
exchange Ions. 

Particle tracking timestep 
constant (dimensionless) 

Velocity diffusion con- 
stant (dimensionless) 

Minimum density to be used 
In calculating barometric 
potential, (m*^) 

Number of Iterative time- 
steps in hydrodynamic 
solution. 


Temperature of Isothermal 
boundary surfaces. (eV) 

Thermal boundary condition 
'ISOT' - Isothermal every- 
where. 

‘SINK* - Isothermal on 
orifice plate. 

'INSU' - Insulating 
orifice plate. 

Radius of ring neutralizer 

Number of potential- 
temperature Iterations. 


OUTPut LUNCHX LUNPDT 

LUNCHX = Logical unit number for charge exhange Ion 
densities and velocities (Default 19). 
LUNPDT = Logical number for current and potential 
printout (Default 20). 

DEST 'NONE' Plot destination 

title 'NASA-S-CUBED Plot title (beginning In 

ION ENGINE Col. 9) 

CODE' 



1 

2 _ 


6 . 

7 

8 . 

9. 

10 . 

11 . 

12 . 

13. 

14. 

15. 

16. 
17. 


NR 24 
HZ 24 
RNAX 15 
ZNAX .30 
CURRENT .OHS 
FLQU .1-* 

SIGMA S.E-19 
TEMP .06 
RADIUS .07 
PROFILE QUADRATIC 
RNEUT .08 
ENERGY 3000. 
THETAI 2. 

NSTEPS SOO 
RHOMIN 1.E10 
THERBC INbU 
END 


Figure 3.1a. 


Input data representative 
of SERT II thruster. 


1 . 

NR 2S 

2. 

NZ 2S 

3. 

RHAX .36 

4. 

ZMAX .72 

5. 

CURRENT 1. 

6- 

FLOU 1.25 

7. 

SIGMA S.E-19 

8. 

TEMP .06 

9. 

RADIUS .15 

iO. 

RNEUT .17 

U, 

ENERGY 1100. 

12. 

NSTEPS SOO 

13 

RHOMIN 1.E10 

1 4 

THERBC SINK 

15. 

END 


Figure 3.1b. Input data representative 
of 30 cm thruster . 



riL''r oPTio'i simnARY 


NR NZ 

RPIAX 

ZHAX 



25 25 

1.0000 1.0000 



RE-VI CURRENT 

BEAM 

RADIUS 

8EAP1 ENERGY 


1.000000 APIPS 

1,1 

DOOO n 

1000.0 F.U 


fAnui 

PROFILE 

SPREAD ANGLE 


200.59 

CAUS 

.0 DEG. 




ION PtASS = 3 

.355—025 KG. 




PEAK CURRENT 

= 5.968 + 018 p1»»t-?3 

<5cc*»i _ 



ION VELOCITY 

= 30893.3 PI /SEC 


SIGPIA 

FLOU 


TEPIP 


6.0-019 

2.000000 

APIPS 

1.000 EU 


ECI1AR 

UDTR 

ETA 

RH0P1IN 

NSTEPS 

10.000 

.200 

.050 

1.00+012 

300 


THFPrtftL BC 


■ISOT' UITH TEHPERATURE 


1.00 EU. 


MeUTO<^I.IZFR RAOIU? = 1.0000 


2 ITERATIONS. 


PRINTED OUTPUT: LOGICAL UNITS NO. 19 AND 20. 


PLOT DESTINATION = NONE 

PLOT title = NASA - S-CUBED ION ENGINE PLASOA PODE 


Figure 3.2a. 


Code Option Summary (Default values). 


I'l DF OPUfU bUnrIArt. 

ilR NZ 

" E 2S 


rtlli-iA ZOAa 

.J6n0 /Vu\ 


BE 'ill rURKCNr BEAPt RADIUS 

:.o0000n AOPS .1500 n 


Btnil ENLPGr 
1100.0 EK 


HAST lAnui 
200 E9 


PROFILE CPPEmO hNCLE 

EAUS ... DEG. 

ION NASS = 3 355-025 I'G. 

PEAK CURRENT = 2.652+020 n*»C-2] 5EC<*t-n 

ION UELOriTY = 32390.7 I1/3EC 


SlfiOA FLOU 

5 '-019 1.250000 AOPS 


TEtIP 
.060 EU 


FCHAR VDTR . ETA 

10 ')uu .200 .050 


RH0I1IN NSIEPS 

1.00+010 500 


THERNAl BC = 'SINK' UITH TEPIPERATURE 


1.00 EU. 


neutralizer RADIUS = .1700 2 ITERATIONS. 


PRINTED OUTPUT: LOGICAL UNITS NO. 19 AND 20. 


PLOT DESTINATION = CALC 

PLOT TITLE = SAPtPLE INPUT DATA 


Figure 3.2b. Code Option Summary (Values similar to 



TABLE 3.2. MAJOR SUBROUTINES OF ION CODE 


Major Subroutines 

(called by main routine) 

Called automatically 

DEFOPT - Sets default options. 

INPUT - Reads user input of problem parameters. 
INPRNT - Summarizes problem parameters. 

Called on user request 


CHANGE [lun]-Change input values by reading from 
logical unit no. lun (default 5). 

CHEX - Hydrodynamic calculation for charge exchange 
ion density and velocity. Results on logical 
unit no. 9. 

RECHEX - Restart of CHEX calculation. 

POTENT - Initialize and iteratively calculate potential 
and temperatures. Results on logical unit 
no. 21. 


NEWPOT - Same as POTENT, but skips initialization and 
begins calculation on new case. 

RESUME - Continues iterative calculation of POTENT 
starting from most recent iterate. 

PLOT - Generate graphical output. 



neutralizer position. The CHANGE subroutine may be used to alter 
problem parameters, and the PLOT subroutine is called for graphical 
output. 


B. INPUT PARAMETERS 

i. Grid Size and Resolution 

The ION code discretizes the ion engine plasma problem on a 
finite element grid of NR by NZ gridpoints. For the charge exchange ion 
calculation each dimension has a maximum of 25. For POTENT the size may 
be increased to 30 by 30. The physical size of the computational space 
is RMAX by ZMAX. It is suggested that ZMAX be approximately double 
RMAX, and RMAX approximately double the beam radius. Resolution should 
be sufficient to model the density falloff at the beam edge. RMAX and 
ZMAX may be reduced between CHEX and POTENT if desired. 

ii. Ion Beam Characteristics 

These variables are self-explanatory, with the exception of 
the beam profile. Two choices are currently available 

GAUSsian; „ 

J(r) = 

QUAD rati c; 

J(r) = Jq 1 - (r/R)^ [r < R] 

where R is the specified beam radius. The GAUSsian profile is 
smoother and gives more reliable convergence, particularly in the CHEX 
calculation. The QUADratic profile is likely to give poor results 
when charge exchange ion production is low. 
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111 . Charge Exchange Ion Production (Physical) 


The rate of charge exchange ion production is proportional 
to 

(a) the cross-section SIGMa, 

(b) the excess of fuel FLOW over beam CURRent; 

(c) the inverse square root of the neutral TEMPerature. 
Typical value for the neutral TEMPerature is 
300“C = .06 eV. 

iv. Charge Exchange Ion Production (Computational) 


These parameters govern the operation of the CHEX 
calculation and, with the exception of NSTEPS, are not recommended for 
change. The "characteristic energy" ECHA, should be roughly the 
barometric potential difference between the beam and the charge 
exchange cloud: 


ECHA e in (, 


beam'^^ch.ex. 


) 


Together with VDTR, this determines the hydrodynamic timestep. The 
velocity diffusion, ETA, is used to promote convergence, and can 
usually be reduced if the timestep is shortened. RHOMIN should 
approximate the minimum density expected in the problem. NSTEPS is 
the number of steps to be run to achieve steady state. (See discussion 
of CHEX routine below). 


V. Thermal and Current Boundary Conditions 

These provide the temperature used in the CHEX calculation, 
as well as the thermal boundary condition and neutralizer position for 
the POTENT calculation. These parameters may be changed prior to 
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executing POTENT or NEWPOT. NITER is the number of Iterations per- 
formed on calling POTENT, NEWPOT, or RESUME. 

v1. Output Destinations 

The printed output from CHEX and POTENT Is rather lengthy. 
For CHEX, the densities and velocities are printed for all node 
points, and for POTENT the densities, currents and potentials. 
Normally, these are sent to scratch files, from which they may be 
edited. To have them actually printed. Input the card 

OUTPUT 6 6 

DEST specifies the site-dependent plot destination on auto- 
matic execution of NASCAP PLOTREAD., and TITLE allows user specifica- 
tion of a plot title, which may be changed during the course of the 
run. 


C. SUBROUTINE CHEX (RECHEX) 

The charge exchange Ion density calculation Is the most 
time consuming part of the ION code, taking 1-1.5 seconds per step at 
S-CUBED. It outputs Immediately the neutral efflux and the resulting 
charge exchange current. Every fifty steps It prints the charge 
exchange Ion density near the beam origin, which gives some Indication 
of convergence. The true test of convergence, however Is the 'TOTAL 
DNDT' printed by the POTENT routine. If this CHEX calculation Is 
properly converged, this number should be within 10 percent of the 
charge exchange current. Otherwise, RECHEX should be used to continue 
the calculation. 

CHEX writes Its Internal Information on logical unit No. 9, 
and prints output on LUNCHX. The printed output lists velocity at 
each node point. The densities printed correspond to positions 
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halfway toward the next R and Z values. Plots of the charge exhange 
ion densities and currents are available after processing by the 
initialization stage of POTENT. 

D. SUBROUTINE POTENT (NEWPOT, RESUME) 

This routine calculates the potentials needed to drive the 
electrons from the neutralizer through the resistive plasma, and the 
temperature profile from the resultant plasma heating. The potential 
boundary condition is that the mean potential is zero. The thermal 
boundary conditions are isothermal at RMAX and ZMAX, and either 
isothermal or insulating at the orifice plate. Alternatively, the 
plasma may be considered isothermal, in which case no heat flow 
calculations are performed. 

After initialization, this routine performs alternate ICCG 
(Incomplete Cholesky Conjugate Gradient) calculations for the potential 
and temperature fields. For each ICCG calculation, convergence is 
indicated by the numerator (RDOTR) being several orders of magnitude 
smaller than the denominator (RDOTRl). For the overall calculation 
convergence requires each denominator being several orders of mag- 
nitude smaller than its initial value. Convergence takes the most 
iterations (~5) for the INSUlating boundary condition. The most 
apparent symptom of incomplete convergence is the appearance of a 
convection cell around the beam edge in the current plots. The RESUme 
command may be used to perform further iterations on an unconverged 
calculation. To perform an additional POTEnt calculation with altered 
thermal boundary condition or neutralizer radius, the NEWPot command 
may be used. 

E. SUBROUTINE PLOT 

This is the user requested subroutine that controls the 
graphical output by ION. It has been designed as an interactive, 
device independent, menu type package but it can also be used in a 
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batch mode of operation. Aside from the required control commands or 
cards, all of the required data is taken from the temporary file 21. 

If that file's data had been saved from a previous run, the user could 
preassign 21., copy the old data to it, PXQT ION, and call PLOT 
immediately to produce or reproduce plots from previous runs. The 
user should be warned however, that if more than one problem is run 
per execution of ION, file 21. will contain only the most recent data 
and the results of the earlier run will have been lost. The output of 
PLOT is coded and stored in file 2. This file should be kept and used 
as input to NASCAP*PLOTREAD which will actually perform the plotting. 
File 2. is not overwritten if multiple runs are made in a single execu- 
tion of ION. Successive calls to subroutine PLOT are automatically 
assigned sequence numbers, and the user may change titles with the 
CHANGE subroutine. 

Upon calling PLOT, the interactive user will see the 
primary menu shown in Figure 3.3a, and will be prompted for the plot 
desired. The appropriate response is the integer associated with the 
desired plot. Following the transmittal of plot request, the menu 
will reappear and another plot may be chosen. This process may be 
repeated indefinitely and is terminated by a blank record. If plot 99 
is chosen, the user is presented the menu shown in Figure 3.3b. These 
are all cross-sectional plots so after responding with the appro- 
priate integer, the user is prompted for the desired Z direction nodal 
index (1,2,3, etc). If the user responds with a negative integer the 
plot will be a constant r cross-section. This process may be repeated 
indefinitely and is terminated by a blank record which will return the 
user to the primary menu. 
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UHICH 


1 

2 

3 

4 

5 

6 

7 

8 
9 
le 
11 
99 
12 
13 

PLOT 


POTENTIAL CONTOURS 
ION DENSITY CONTOURS 
ION DENSITY (LOO 
ION CURRENT HECTORS 
ION CURRENT (LOO 
TOTAL CURRENT VECTORS 
TOTAL CURRENT (LOO 
TOTAL CURRENTXRADIUS 
ELECTRON CURRENT VECTORS 
ELECTRON CURRENT (LOO 
ELECTRON TENPERATURE 
RADIUS PLOT 

CHARGE EXCHANGE ION DENSITY CONTOURS 
CHARGE EXCH«(GE ION CURRENT CONTOURS 
>99 


Figure 3.3a. 


CHOSE CROSS- 
1 
2 
3 

UHICH PLOT > 


SECTIONAL PLOTS 
POTENTIALS 
ION DENSITY 
ELECTRON TENPERATURE 


Figure 3.3b. 


Figure 3.3. Subroutine PLOT menus 



4. TEST CASES 


In this Chapter we present four test cases for parameters 
representing the SERT II thruster, and two calculations for a 30 cm 
diameter thruster. These calculations were performed prior to 
collecting and refining the code to its present form, so we omit any 
detailed output. Also, the plots presented were made using the 
DISSPLA plot package and differ somewhat from the NASCAP-like plots 
produced by the final version of the code. Nonetheless, we believe 
these results represent well the type of calculation performed by the 
ION code. 


A. SERT II CALCUUTIONS 

The SERT II thruster was modeled for the ION code as having 
a 7 cm radius ion beam. The beam had consisted of 0.085 amperes of 
3 keV Hg ions. The beam profile was assumed to be 

J (r,z) = Jg (1 - (r/.07)2) 

An additional 0.055 ampere equivalent flow of 0.06 eV neutrals 
resulted in 0.9 mA of charge exchange ions, which expanded in 
barometric potentials characterized by a 2 eV electron temperature. 

The thermal boundary conditions at RMAX = 0.20 and ZMAX = 0.40 were 
set to 2 eV, while the orifice plate could be either insulating or a 2 
eV heat sink. The neutralizer radius was either 8 cm or 10 cm. The 
results for maximum potential and temperature are shown in Table 4.1. 
It is apparent that the orifice plate boundary condition has 
substantial effect, and also, that moving the neutralizer closer to 
the beam has some effect in reducing joule heating of the plasma. 
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TABLE 4.1. SERT II RESULTS 


Neutral izer 
Position 

Thermal Boundary 
Condition 2 eV 

Insulating Thermal 
Boundary 

0.08 

& = 2.3 eV 

max 

i = 13.8 V 

max 

0 = 4.1 eV 

max 

1 = 24.2 V 
max 

0.10 

e = 2.4 eV 

max 

«'max = W-' ' 

0 = 4.5 eV 

max 

6 = 24.2 V 

max 


The plasma density for this model is shown in Figure 4.1. 

15 -3 

The density varies from over 10 m at the beam center to 
~10^^ m’^ at the upper right of the plot. The sharp drop at the 
beam edge is apparent. It is also clear from the figure that near the 
orifice plate the beam plasma density and charge exchange plasma 
density are comparable. 

Figures 4. 2-4. 5 show the results of these runs for 
potential, electron temperature, and total current. In all cases, it 
appears that the hottest temperatures occur near the neutralizer. The 
resemblance between the potential contours and logarithmic density 
contours indicate the extent to which the barometric law is valid. 

The total current plots indicate beam neutralization within about one 
thruster diameter of the orifice; downstream current patterns result 
from incomplete convergence. There is some indication of a shorter 
neutralization length for the 2 eV orifice plate than for the 
insulating orifice plate. 
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B. 30 CM THRUSTER CALCULATIONS 


Two further calculations were run comparing insulating and 
heat sink boundary conditions for a 30 cm diameter thruster. The 
thruster emitted one ampere of 1100 eV Hg ions, and produced 20 mA of 
charge exchange ions. Thus the current densities and plasma densities 
are both appreciably higher than for the SERT II cases. The ambient 
electron temperature was set to 1 eV, and the neutralizer radius to 
17 cm. 


The beam profile used was a Gaussian: 

J{r) = exp C-3(r/rgg)^] 

where the nominal beam radius, rgg, contains 95 percent of the beam 
current. The resulting charge exchange ion density (Figure 4.6b) and 
total density (Figure 4.6a) is much smoother than in the SERT II case, 
where a quadratic beam profile was assumed. 

Despite the lower electron temperature boundary conditions, 
the results (Figures 4. 7-4.8) show temperatures and potentials similar 
to, or even greater than, the comparable SERT II cases. For the 
insulating case, we see the temperature maximum at the beam center 
rather than near the neutralizer. Otherwise, the qualitative 
conclusions are similar to those of the previous cases. 
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(a) 


(b) 


(c) 


Figure 4.3. Temperature (a), current density (b), and potentials (c), for SERT II case: 2 eV orifice 

plate; 10 cm neutralizer radius. 
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5. SAMPLE RUN 


In this Chapter we present a sample run of the ION CODE. The 
run was performed in batch mode on S-CUBED's UN I VAC 1100/81, using 
ASCII FORTRAN Level 9R1 with full optimization, and the plots were 
done on the Gould electrostatic plotter. The run conditions represent 
variants on the SERT II thruster. 

Figure 5.1 shows the runstreara for the sample run. The initial 
portion of the runstream involves breakpointing to a print file, 
printing a header and printing the runstream. After invoking the ION 
code C@XQT BIGION.BIGION] the initial input is OADDed from data 
element BIGION.SERT. The remaining input is explicitly shown. At the 
end of the runstream, restart files 9 and 21 are copied from the 
temporary files created by ION into permanent files. (Alternatively, 
ION could have run directly on permanent files by 


0ASG,PU SAVE 9. 

0ASG,PU SAVE 21. 
OUSE 9, SAVE 9 

OUSE 21,SAVE21.) 


Figure 5.2 shows the first page of output from the ION code. 

ION issues a welcome message and requests input. The input cards 
(resulting from the OADD statement) are echoed. Then a complete 
option summary is printed. Next the charge exchange ion calculation 
(CHEX) is invoked. The neutral efflux and total charge exhange cur- 
rent are printed. As this hydrodynamic calculation is rather lengthy, 
a message is printed every fifty steps. On invoking POTENT, we notice 
that the 'TOTAL DNDT' does not match the charge exchange production, 
indicating incomplete convergence in CHEX. The four ICCG messages 
constitute two iterations of potential and temperature calculations. 
(The GETCH warning may be ignored if the subsequent fraction indicates 
convergence.) Convergence is achieved when, in the final pair of 
ICCG's the numerators are small compared with the denominators, and 
(more important) the denominators are small compared with their 
initial values. 
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Figure 5.3 shows the next section of output. The PLOT 
subroutine was requested, and we were presented with a plot menu. We 
requested plot 1 (potentials). The plot menu repeats (not shown) 
after each plot. We continued by requesting plots 3, 12, and 13, 
terminating with a blank card. The plots, together with title frame, 
are shown in Figures 5. 4-5. 8. 

The printout continues in Figure 5.9. We prepare to restart the 
CHEX calculation by changing the number of steps to 50. While about 
it, we route the plots to the electrostatic plotter (S-CUBED feature). 
RECHEX is called to continue the hydrodynamic charge exchange calcula- 
tion, and DNDT, as calculated by POTENT, is significantly better. We 
then request again the same plots (Figures 5.10-5.14; plot menus 
omitted. ) 

In Figure 5.15 we repeat RECHEX and POTENT, now having rather 
good convergence. The third sequence of plots appears in Figures 
5.16-5.20. In Figure 5.21 we RESUME the potential calculation where 
we left off in 5.15. This time we plot potentials, temperatures, and 
currents (Figures 5.22-5.25). The convection cell around the beam 
edge (Figure 5.25) is indicative of incomplete convergence. 

In Figure 5.26 we change the thermal boundary condition to 
ISOThermal. NEWPOT is invoked for a potential calculation using the 
same ion densities as previously. We plot potentials and currents 
(Figures 5.27-5.29). In Figure 5.30 we repeat the calculations for a 
1 eV orifice plate. Here we see good convergence in four iterations. 
Plots of potential, temperature and current are requested (Figures 
5.31-5.34). Finally (Figure 5.35) an END card exits us from the ION 
code, we execute the DISSPLA plot interface routine, and complete the 
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ac 

eOF *8C 


HJ,1.11idd-UU.riANULLL'-nv^Uf>UI. 

SOruCTCfC &1G10NTCST. 
aASbfUP 8IGI0NTeST* 

9FPCC BIGI0NTCS1* 
ftASG.A alGIONTtST, 

SBRkPT PPIhM/BIGlONTtST 
aalb-NAMC.P 

lOM cooc run 

eox 24 

SnCGiN vbid •Mfd8f4y4 
(»£0»R RuN$TAPT*eiGlON 
LNP* 

aXOT blGlON»BIGION 
8A0D BIGIOK.SCPT 
CHEX 
POTENT 
PLOT 
1 
3 

12 

13 

CHANCE 

NSTEPS 50 

DE5T ELECTROSTATIC 

END 

RECHEA 

POTENT 

PLOT 

1 

3 

12 

13 

PECHEX 

POTENT 

PLOT 

I 
3 

12 

13 

RESUME 

PLOT 

1 

II 
8 

CHANGE 
THER6C ISOT 
END 

KEaPOT 

PLOT 

1 

6 

CHANGE 
THERBC SINK 
NITER 4 
END 

NEaPOr 

PLOT 

1 

11 

B 

ENO 

SPMUiEL 

aASbfPU SAVE9* 

8ASG,PU SAVE21. 
aCOPY 9. ,SAVE9 • 
aCOPY 21.,SAVE21. 
aaiG-NAME.p 
ION CODE Aun 
ENO 

BOX 24 

5BRKPT PRINTS 
SFREEM. 

a$YH,U aiGIONTEST. 

GFIN 


NO CORRECTIONS APPLIEOt 


iXOT BIGION.PIGION 


Figure 5.1. 


Runstream for sample run, as printed by lines 
stream, followed by @XQT command invoking ION 


11-12 of run- 
code. 
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WELCOHE TO ION COMPUTER CODE 
0»TEiC8118l TIMErilJJST 

<PtE*SE ENTER OR 3«00 INPUT P»R*METERS> 

NR 2<« 

N2 2* 

RN«X .15 
2NAX .10 
CURRENT .085 
FLOa .!<■ 

SIGH* 5.E-19 
TEMP .05 
RADIUS .07 
RNEUT .085 
ENERGY 3C00. 

NSTEPS JDO 
RHOHIN l.EIO 
THERBC INSU 
END 

CODE OPTION SUMMARY 


NR N^ 

RHAX 

ZMAX 


21 

.1500 .3000 


BtAH CURRENT 

BEAM 

RADIUS 

BEAM energy 

•C8S000 AHRS 

• 1 

0700 M 

SOOQ.O EV 

MASS (AMU) 

PROFILE 

SPREAD ANGLE 

200.59 

GAUS 

«C DEG. 



ION HASS : 

3.355-025 KG. 



PEAR CURRENT 

= 1.035*020 MA* 



ION velocity 

= 53R91.5 M/SE' 

SIGMA 

FLOW 


TEMP 

5.0-019 

.190000 

AMPS 

•060 EV 

£CHAR 

VO TR 

eta 

RHOMIN 

10.003 

.200 

• 050 

1.00*010 

THERMAL BC = * 

INSU* WITH 

TEMPERATURE 

1.00 EV. 


■21 SEC**«-l) 


NSTEPS 

JQO 


NEUTRALI2ER RADIUS : .0850 

PRINTED OUTPUT: LOGICAL UNITS 


PLOT DESTINATION 
PLOT TITLE r 


: NONE 

NAS* - S-CU8E0 


2 iterations. 

NO. 19 AND 20. 

ION ENGINE PLASM* CODE 


••♦time LEFT : 1792 SEC0N0S»*9 


•«*»**CHEX 

NEUTRAL EFFLUX : 


.0550 EOUIVALENT AMPS. 


total 

GENERATION 

RATE : .000773 

AMPS 



2.906*014 

CHEX - 

50 

STEPS 

completed. 

RHOI 

i .OOSy 

• 007) 

z 

CHEX - 

100 

STEPS 

COMPLETED. 

RHO< 

1 .OOSy 

• 007) 

z 

4.43e'»Ql4 

CHEX - 

150 

STEPS 

completed • 

RHOI 

[ .COS, 

.007) 

z 

9«82S*01<4 

CHEX - 

200 

STEPS 

COMPLETED. 

RHO( 

• 005, 

• G07) 

z 

4 .764*014 

CHEX - 

250 

STEPS 

COMPLETED. 

RHOI 

[ .COSy 

• C07) 

z 

4 .665*014 

CHEX - 

300 

STEPS 

COMPLETED. 

RHOI 

: .005, 

.007) 

- 

4.565*014 

•••TIME 

LEFT 

: 1391 SEC0N0S*pp 






»OTENT 


CURRENT ENTERING FROM LEFT : 3.59-002 

CURRENT EXITING AT RIGHT : 8.59-002 

CURRENT EXITING AT RMAX : 5.89-009 

TOTAL DNDT : -5.8«-00» 

EMITTER AT NODE 19. 

GETRHS — CURZO.ELJ : 8. 59-002 2.96*001 


••*HARNING*«9 GETCH — - CH( 
ICCG RDOTR/ROOTRl 


5977) I -1.00*000 
2.32-005/ 7.75-001 


ICCG RDOTR/R007R1 r 3.18-019/ 1.35-001 

•••hARNlNG*** GETCH CH( 59771 = -1.00*000 


niter: 31 
NITER: 25 


ICCG -— ROOTR/ROOTRI : 
ICCG RDOTR/ROOTRI : 


2.16-003/ 5.79-001 

7.92-015/ 9.92-003 


niter: 31 

NITER: 29 


Figure 5.2. Output resulting from lines 13-16 of runstream. 
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1 POTENTIAL CONTOUPS 

2 ION density contours 

I ION density ilogi 

P ION CURRENT VECTORS 

5 ION CURRENT ILOG) 

6 TOTAL CURRENT VECTORS 

7 total CURRENT ILOGT 

3 TOTAL CUHRENTVRAOIUS 

? ELECTRON CURRENT VECTORS 

1C ELECTRON CURRENT (LOG) 

II ELECTRON TEMPERATURE 

99 RADIUS PLOT 

1? CHARGE exchange ION DENSITY CONTOURS 

1! CHARGE EXCHANGE ION CURRENT CONTOURS 

WHICH PLOT XV 


gure 5.3. Plot menu (line 17 of runstream) 



NASA - 3-CUBED ION ENGINE PLASMA CODE 


DATE-08/ 11/81 TIME-ll. ■13:08 


I 



Figure 5.4. Plot title page. 
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SIXU-Z 


PLOT SEQUENCE » 


1 


NflSfl - S-CUBED ION ENGINE PLflSNfl CODE 


POTENTIRLS (VOLTS) 



R-flXIS 


MIN - 
MRX - 


CONTOUR VALUES 
-J. 0000+000 


-3.4347+000 

1.4530+001 


-3.0000+000 
-2.0000+000 
-1.0000+000 
-7.4506-008 
1 0000+000 
2.0000+000 
3.0000+000 
4 0000+000 
5.0000+000 
6.0000+000 
7 0000+000 
3 0000+000 
9 0000+000 
1.0000+001 
1. 1000+001 
1.2000+001 
1.3000+001 
1.4000+001 


Figure 5.5, Potential contour plot requested in line 18 
of runstream. 
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Z-fiXIS 


PLOT SEQUENCE » 


1 


NRSfl - S-CUBED ION ENGINE PLflSm COOE 


LOG ION DENSITY 



R-fIXIS 


MIN - 1.2000+001 

MBX - 1.5376+001 


CONTOUR VALUES 
1.2000+001 
1.2200+001 
1.2400+001 
1 2600+001 
1.2800+001 
1.3GOO+noi 
1.3200+001 
1.3400+001 
1.3600+001 
1.3800+001 
1.4000+001 
1.4200+001 
1.4400+001 
1.4600+001 
1,4800+001 
1.5000+001 
1 5200+001 


Figure 5.6. Logarithmic ion density contour plot requested 
in line 19 of runstream. 
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Z-RXIS 


PLOT 3EQIJENCE * - 1 NflSR - 3-CUBED ION ENGINE PLASMA COOE 

CHEX ION DENSITIES 



R-AXIS 


CONTOUR VALUES 
1.0000-037 
5.0000+013 

MIN - 9.-I667+0U 1.0000+014 

I. 5000+01 4 
2.0000+014 

MAX - 4.6942+014 2.5000+014 

3.0000+014 
3.5000+014 
4.0000+014 
4.5000+014 


Figure 5,7. Charge exchange ion density contour plot requested 
in line 20 of runstream. 
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»»IIME LEFT 


13>lo SECONDS*** 


«**«*«CHANCE 

NStEPS 50 
OEST ELECTROSTATIC 
***FBREAK**» CROUP NO. 
OEST ELECTROSTATIC 
ENO 

CODE OPTION SUMNAPY 


2 IN THIS CARO IHAGE - TOO LONG 


NR 

2N 


N2 

2<t 


RMAX 

.1500 


2HAX 

.3000 


beam current 
.035000 AHPS 

HASS lAHU) 
200.59 


beam RADIUS 
.uTOC H 


BEAM ENERGY 
3000.0 EV 


PROFILE SPREAD ANGLE 

GAU5 .0 OEG. 

ION HASS : 3.355-025 KG. 

peak current : 1.035*020 H**t-21 SEC**t-ll 


ION velocity : 53N91.5 H/SEC 


FLOW 

.INOOOO 


lEHP 
.060 EV 


VDTR 

.200 


SIGMA 
5.0-019 

ECHAR 

10.000 

THERhAl BC : 'INSU* WITH TEHPERATURE 
NEUTRALI2ER RADIUS = .0BS3 

PRINTED OUTPUT, logical UNITS NO. 19 

PLOT TITlE^:^^°** ' NaIa - S-CUBED ION ENGINE PLASMA CODE 
***TIME LEFT : 13Ao SECONDS*** 


ETA RHOMIN 

.050 1.00*010 

1.00 EV. 

2 ITERATIONS. 

ANO 20. 


NSTEPS 

50 


neutral efflux : .0550 equivalent AMPS. 

TOTAL GENERATION RATE = .000773 AMPS 

CHEX - SO STEPS COMPLETED. RHOt .005, 


***T1ME LEFT 


1275 SECONDS*** 


.007) : 


4 .n9S*014 


******POTENT 

CURRENT ENTERING FROM LEFT = 8.54-002 

CURRENT EXITING AT RIGHT : 8.54-002 

CURRENT EXITING AT RHAX : 6.92-004 

TOTAL ONOT = -6.93-004 

EMITTER AT NODE 14. 

GcTRhS — CufiZO.ELJ X 8. 54-002 2.46*001 

•**WARNING*** GETCH CH( 5477) : -1.00*000 

ICCG RDOTR/ROOTRl r 2.30-002/ 7.75-001 

ICCG ROOTR/ROOTRl : 1.34-013/ 1.41-001 

*.*.ARNING*** GETCH CH( 5477) = -1.00*000 


niter: 31 
niter: 24 


2.16-006/ 7.14-001 niter: 31 
RDOTR/ROOTRl I 7.73-015/ 9.73-003 NITER: 23 
: 1233 seconds*** 


ICCG RDOTR/ROOTRl : 

ICCG - — 

***T1ME LEFT 

******PLOI 


Figure 5.9. Output resulting from lines 23-29 of runstream, 



NflSfl - S-CUBED ION ENGINE PLASMA CODE 


OATE-08/11/0] TIME-U: 45:31 



Figure 5.10. Plot title page (line 29 of runstream). 
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sixy-z 


PLOT SEQUENCE >» 


2 


NASA - S-CUBED ION ENGINE PLASMA CODE 


POTENTIALS (VCLTS) 



MIN 

MAX 


CONTOUR VALUES 


-4.3982+000 

1.3340+001 


-5.0000+000 
-4 0000+000 
-3.0000+000 
-2.0000+000 
-1.0000+000 
-1.3411-007 
1.0000+000 
2.0000+000 
3.0000+000 
4.0000+000 
5.0000+000 
6.0000+000 
7.0000+000 
8.0000+000 
9.0000+000 
1 0000+001 
1.1000+001 
1.2000+001 
1.3000+001 


P-AXIS 


Figure 5,11. Potential contour plot (line 30). 
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PLOT SEQUENCE » - 2 NflSfl - 3-CUBEO ION ENGINE PLPISm CODE 

LOG ION DENSm 



MIN 

I1RX 


R-flXIS 


1 2015+001 
1.5374+001 


CONTOUR VRLUES 
1.2000+001 
1 2200+001 
1.2400+001 
1.2600+001 
1.2800+001 
1.3000+001 
1.3200+001 
1.3400+001 
1.3600+001 
1.3800+001 
1.4000+001 
1.4200+001 
1.4400+001 
1.4600+001 
1.4800+001 
1 5000+001 
1.5200+001 


Figure 5.12. Ion density logarithmic contour plot (line 31) 



PLOT SEQUENCE « - 2 


NflSn - S-CUBED ION ENGINE PLflSHfl CODE 


CHEX ION DENSITIES 



1.0339+012 


4.6497+014 


CONTOUR VALUES 
1.0000-037 
5.0000+013 
1.0000+014 
1.5000+014 
2.0000+014 
2.5000+014 
3.0000+014 
3.5000+014 
4.0000+014 
4.5000+014 


.000 .025 .050 .075 .100 125 

R-RXIS 


Figure 5.13. Charge exchange ion density contours (line 32). 





•**TIHE LEFT 


1231 SECONDS*** 


«**«**RE CHEX 
NEUIRiL EFFLUX : 
total OENERATION 
CriEX - 5D STEPS 

***TIME LEFT : 11 


•□55D equivalent amps. 

RATE = .000773 AMPS 

completed. RH0( .005, 

60 SECONDS*** 


.DOT) = 


V .NSS^OIN 


«*****P0 TENT 

CURRENT ENTERING FROM LEFT : 8. 
CURRENT EXITING AT RIGHT : a. 
CURRENT EXITING AT RMAX 


total DNOT 


= - 1 . 


EM 

ITTER 

AT 

NODE 


14. 



GE 

TRHS - 

• 

CURZO. 

EL 

J = 6 

.5R-C02 

«« 

vmARNI 

NG 

*** G& 

TC 

H — 

CH( 

5477 


ICCS 

- 

— HOO 

TR 

/ROOTR 

1 = 

!• 


ICCG 

• 

— ROO 

TR 

/ROOTR 

1 = 

1. 


•kARNI 

NG 

*** Gc 

TC 

H — - 

Ch( 

5477 


ICCG 

- 

— ROO 

TR 

/ROOTR 

1 = 

1. 


ICCG 

- 

— ROO 

TR 

/ROOTR 

1 = 

4. 


• TIME 

LC 

FT = 

11 

18 SECOMOS««» 


5V-002 

SM-0D2 

N J-QQIt 

rj-COR 

.>46«C01 
: -1.00»ODO 

-OOR/ 7.77-001 

-013/ 1.33-001 

: - 1 . 00*000 

-OCR/ 6.21-001 

-015/ 8.36-003 


NITER: 31 
NITER: ZR 

niter: 31 
niter: 2R 


PLOT 


Figure 5.15. 


Output resulting from lines 35-37 of runstream. 
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NfISfl - S-CUBED ION ENGINE PLflSMB CODE 


DflTE-08/11/81 TIME-ll:-)8:22 


Figure 5.16. Plot title page (line 37 of runstream). 
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Z-fiXIS 


PLOT SEQUENCE » 


3 


NRSR - S-CUBED ION ENGINE PLfiSMfl CODE 


POTENTIRLS (VOLTS) 



R-flXIS 


CONTOUR VALUES 
-5.0000+000 
-■I. 0000+000 

MIN - -4.3779+000 “3.0000+000 

-2.0000+000 
-1. 0000+000 

MAX - 1.2881+001 -1.3411-007 

1. 0000+000 
2.0000*000 
3.0000+000 
4.0000*000 
5.0000+000 
6.0000+000 
7.0000+000 
a. 0000+000 

9.0000+000 

1.0000+001 

1.1000+001 

1.2000+001 


Figure 5.17. Potential contour plot (line 38). 
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PLOT SEQUENCE « - 3 


NASA - 3-CUBED ION ENGINE PLASMA CODE 


LOG ION QENSm 



R-AXIS 


CONTOUR VALUES 
1.2000+001 
1.2200+001 

MIN - 1.2021+001 

1.2600+001 

1.2800+001 

MAX - 1.5373+001 1.3000+001 

1.3200+001 
1 3400+001 
1.3600+001 
1.3800+001 
1.4000+001 
1.4200+001 
1 4400+001 
1.4600+001 
1.4800+001 
1.5000+001 
1.5200+001 


Figure 5.18. Ion density logarithmic contours CT^ine 39) 
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Figure 5.19. Charge exchange ion density contours (line 40). 
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Figure 5.21. Output resulting from lines 43-44 of runstream. 
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Figure 5.23. Potential contour plots (line 45). 
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Figure 5.24. Electron temperature contours (line 46). 
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Figure 5.26. 


Output resulting 


from lines 49-53 of runstream. 
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Figure 5.27. Plot title page (line 53). 
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Figure 5.28. Potential contours (line 54). 
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Figure 5.30. 


Output resulting 


from lines 57-62 of runstream. 
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Figure 5.32. Potential contour plot (line 63). 


62 



Z-FlXIS 


PLOT SEQUENCE b - 6 


NASA - S-CUBED ION ENGINE PLfiSm CODE 


ELECTRON TENPE3IflTURE3 



.000 .025 .050 .075 .100 .125 .150 

R-flXlS 


MIN - 
MRX - 


1.0000+000 

1.3505+000 


CONTOUR VALUES 
9.8000-001 
1.0000+000 
1.0200+000 
1.0400+000 
1.0600+000 
1.0800+000 
1.1000+000 
1.1200+000 
1. MOO+000 
1.1600+000 
1.1800+000 
1.2000+000 
1.2200+000 
1.2400+000 
1.2600+000 
1.2800+000 
1.3000+000 
1.3200+000 
1.3400+000 


Figure 5.33. Electron temperature contour plot (line 64). 
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Figure 5.35. Output resulting from lines 67-72 of runstream, 
including execution of graphics interface to 
DISSPLA plot package. 
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Abstract } 

The purpose of the present study is to 
determine the capability of a fluid model 
of electron transport to explain observed 
properties of ion thruster generated plas- 
mas. Calculations reported here show that 
when the effective collision frequency in 
such a nodal is of the order of the elec- 
tron plasma frequency, the resulting elec- 
tric potential variations and electron 
temperatures are in qualitative agreement 
with values measured in the plasma gene- 
rated by the SERT II thruster. Both theory 
and probe measurements made in flight and 
ground tests indicate substantial depar- 
tures from the barometric law and strong 
variations of plasma potential across the 
beam boundary. 

Nomenclature 
3 nagnetic field 

C collision operator in Boltzmann- 

Vlasov equation 
S _ _ electric field 

f^drdv number of particles of type a in 
phase_^spaoe volume element dr dv 
at r,v 

? total energy flux 

F^ = q^(E vxB/c) 

I unit tensor 

3 net current density 

k Boltzmann's constant 

mass of paruicle of species a 

m electron mass _ _ _ 

n electron density =• /f(r,v,t)dv 

p = 1/3 n /f v'^dv = scalar electron 

pressure 

? m J"f v'v'dv 3 pi + pressure 

tensor , 

Q = 1/2 m/v'^^C dv 

charge on particle of species a 

q magnitude of electron charge 

3 heat flux, Eq. (7) 

? position vector of a particle 

R 3 /mv C dv 

T electron temperature 

y velocity of a particle 

v' 3 V - V 

V mean or_drift velocity » 

1/n Jf(r,v,t)dv 

' mean free path for pair collisions 

between electrons 

0 3 icT 

" plasma resistivity 

< thermal conductivity of plasma 

. effective collision frequency 

electron-ion collision frequency 

T 3 nm TV - '"V'^>/3 I = stress tenson 


electric potential 
electron plasma frequency 

1 . Introduction 

The purpose of the present study is to 
determine the capability of a fluid model 
of electron transport to explain observed 
properties of ion thruster generated plas- 
mas. Calculations reported nere show that 
when the effective collision frequency in 
such a model is of the order of the elec- 
tron plasma frequency, the resulting elec- 
tric potential variations and electron 
temperatures are in qualitative agreement 
with values measured in the plasma gener- 
ated by the SERT II thruster. Probe mea- 
surements made in SERT II flignt and 
ground test experiments (2 ) indicate sub- 
stauitial departures from tne carometric 
law<3-7) show strong variations of 

plasma potential across the beam boundary. 

We propose to explain t.he plasma proper- 
ties observed in the aforementioned experi- 
ments in terms of anomalous resistance of 
the plasma to t.he flow of electron current. 
The calculations are based on fluid equa- 
tions expressing conservation of charge, 
momentum, and energy. We adopt the clas- 
sical (ignoring thermoelectric effects) 
form of. the equations of electron trans- 
port, but permit reduced values of the 
transport coefficients. Predicted space 
dependent potentials, electron tempera- 
tures and current densities agree qualita- 
tively with experimental results. 

While the plasma is not collision domi- 
nated, randomization of electron velocities 
may still occur tnrough en.nanced levels of 
fluctuating fields, such as those initiated 
by streaming instabilities. Such fields 
are probably effective in coupling neutral- 
izer electrons into the bulk plasma and in 
equalizing the mean drift of electrons wit.n 
ions in the thruster beam. Such mec.nanisms 
are often approximated by introducing an 
effective collision frequency, .. 

Within a meter or so of the ion thruster, 
the electron densities are in the range 

10^ < n < 10^^ cm"^ 

and their velocity distribution is charac- 
terized by temperature 5 between about one 
and ten electron volts. The Debye length 



IS typically small compared to distances L 
over which there is a substantial variation 
of macroscopic plasma properties such as 
density, potential, and temperatures. On 
the other hand, the mean free path for pair 
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collisions 

X a 10“^ 3^^^/n ca, E < 3 

c *• 


for electrons of energy S (eV) is typically 
long compared to L, so that as previously 
asserted the behavior of the plasma is con- 
trolled by collective rather than colli- 
sional effects. Since Xq << L, the plasma 
IS quasineutral, departures from neutrality 
amounting roughly to 


an/n 



10 


-4 


the space around the vehicle is strongly 
shielded from surface potentials. This is 
in contrast to the situation that prevails 
in charging of spacecraft in geosynchronous 
orbit where affects of space charge are 
entirely negligible and potentials are 
determined as solutions of Laplace's equa- 
tion. 


Although collisionless, thruster-gene- 
rated plasmas exhibit Tiacroscopic behavior 
similar in many respects to that of a col- 
lisional plasma. Such behavior is perhaps 
not totally unexpected in view of the fact 
that in both non-equilibrium and equili- 
brium plasmas electrons are scattered by 
fluctuating electric fields. A primary 
difference between the equilibrium and non- 
equilibrium cases is in the magnitude of 
the fluctuating fields. 


Several investigators have measured 
properties of thruster generated plas- 
mas. In the experiments of Ogawa, 

at al. , on cesium ion oeams neutralized 
by electrons from a hot wire, measurements 
were made of the density, potential, and 
electron temperature in the beam plasma. 
The potential difference between the 
neutralizer wire and the plasma could be 
varied by changing the position of the 
wire, the large potential differences 
(electron iroection voltages) occurring 
when the wire was completely withdrawn 
from the beam plasma. An important result 
of the Ogawa exoeriments was that over a 
wide range of conditions electron density 
n and plasma potential 5 were well corre- 
lated by the barometric law 


n(r) = const exp (qo (r) /kT) 


( 1 ) 



Since the barometric law is a thermal 
equilibrium concept, it can be completely 
valid only if the plasma is isothermal. 

The plasma is only approximately isotherm- 
al, noticeable deviation occurring as one 
proceeds from the beam axis beyond the 
beam edge into the plasma formed by am- 
bient and charge exc.hange ions. Kaufman 
observes an electron temperature in the 
charge exchange plasma only about half 
that in the beam. Ogawa and 
Sellen(2' obtained measurable temperature 


variations in the beam plasma over several 
tens of centimeters in the downstream direc- 
tion from the accelerator grid. The largest 
deviations from the barometric law were 
observed for large infection potentials 
(^10 volts) . Probe traces in such cases 
also indicated departures of the electron 
spectrum from a Maxwellian shape. 

Probe measurements of the plasma poten- 
tial in the thruster beam were made in 
SEHT II flight and ground test experiments. 
The measurements show strong variation of 
plasma potential across the beam boundary 
about 20 cm downstream from the thruster 
grids. Such results are difficult to ex- 
plain on the basis of a barometric law 
relationship unless the electron tempera- 
ture or density variation from beam center 
to beam edge is much nigher than might be 
expected from other measurements made in 
similar configurations, we anticipate, 
however, that such is not the case and, in- 
stead, that the observed behavior should 
be explained in terms of the anomolous 
resistivity of the thruster generated plas- 
ma to the flow of electron current. Thus, 
the primary objective of the following 
sections of the report is to determine the 
capability of simple transport models to 
explain, at least qualitatively, the ex- 
perimental results. 

The next section summarizes the ki.netic 
equation for the electron distribution and 
the first few moment equations expressing 
conservation of charge, momentum, and 
energy. In Section 3 we state the approxi- 
mations leading to the transport aquations 
which we eventually solve. The method of 
solution and the results of calculations 
will be the subject of Sections 4 and 3 
respectively. The final section. Section 6, 
summarizes t.he conclusions of this study. 

2. E.xact Scuations for Electron Gas 


In principle, a kinetic approac.h based 
on the Vlasov-3oltzmann equation fully 
describes the spacecraft ge.nerated plasma. 
The complexity of such an approach, How- 
ever, makes it impractical as a basis for 
conducting multidimensional calculations 
of plasma behavior. Besides, except near 
sources and collecting surfaces, where the 
distribution function may change markedly, 
one should be able to adequately describe 
the plasma in terms of certain average 
properties of the distribution, such as 
temperature, density, and particle and heat 
fluxes. Below, the exact equations de- 
scribing the plasma are given in order that 
the reader may be aware of the effects 
neglected in arriving at the approximate 
equations that are subsequently solved. 

Quite generally the state of the plasma 
can be specified by the distribution func- 
tion fa(r,v,t) that characterize each parti- 
cle component a, where £^(c,v,t)<ir dv 
represents the number of particles of 
species a_i.n_the six dimensional volume 
element dr dv about the position r,v in 
phase space. The kinetic equations whic.h 
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describe the distribution are 




7 f 
V a 


( 2 ) 


For particles of aass and charge in 
an electric field £ and a magnetic field 
3 the "smoothed" force on a particle is 
Pj. The effects of collisions between 
particles is taken account of by the col- 
lision term denoted here by C^. Here we 
attempt to describe the plasma in terms of 
Its density a, mean velocity V, and cer- 
tain nigher velocity moments. For conven- 
ience, we have omitted the particle 
species subscript a. The first three 
moments of the kinetic equation yield 
conservation equations for particles, 
momentum and energy, as summarized be- 
low: O) 


Conservation of Particles 

^ + T-nV = 0 (3) 

Conservation of Momentum 

1^ (mnV) + 7-mnW + 7*f - qn/S + 

(4) 

Conservation of Energy 

^ , qu2.v 

R*V + Q (5) 


electrons as a whole and that their distri- 
bution varies slowly in space. 

3. Approximations for Electron Gas 

Consider that the plasma is in a steady 
state and that quasi-neutrality pertains 
throughout the bulk plasma (that is, away 
from electrodes and collecting surfaces) . 
The electrons and ions each satisfy the 
particle continuity equation 

7*nj_V, =0 (1 =• +,-) (9) 

with m ■ n_ » n. The momentum equation 
simplifies considerably if the electron 
drift velocity V is small compared to the 
random velocity <v'^> 2/2 and if the veloc- 
ity distribution is nearly isotropic. 

Then, in the absence of magnetic fields, 

7p + en£ =• R (10) 

where R represents the collisional drag be- 
tween ions and electrons. In a classical 
plasma dominated oy collisions, R is com- 
posed of a part rela- 

tive motion u = Vg - between electrons 
and ions, leading’to plasma resistivity, 
and to a thermal part proportional to the 
gradient of electron temperature, which is 
frequently neglected. In this approxima- 
tion, equation (4) becomes 

7p t enE =• ine 3 (11) 

where 3 is the net current density and the 
plasma resistivity ^ is related to the 
electron-ion collision frequency vgj_ by 


where 




LS the total energy flux. 


c » nm 




V 


T • V * q 
( 6 ) 


(7) 


IS the heat flux, 

r -2 

Q 


(3) 


and < > denotes an average over t.he dis- 

tribution f. 


So far, the equations are quite general 
and involve no assumption that the gas is 
collision dominated or retains a Maxwellian 
spectrum of velocities. Separate conserva- 
tion equations may be written not only for 
different particle species, but also for 
different groups or particles of the same 
charge and mass. Primary electrons, for 
example, with significant streauning ener- 
gies could be treated as distinct iron t.he 
mam electron population which is taken to 
have a .Maxwellian distribution of veloci- 
ties. For the present, however, and until 
experimental or theoretical considerations 
dictate otherwise, we shall consider 


n 


-1 


p 

4v 


ei 


( 12 ) 


If the plasma is non-resistive and iso- 
thermal, equation (11) yields the baro- 
metric law, aquation (1). In this sense, 
aquation (11), or more generally the com- 
plete electron momentum equation, may be 
regarded as the generalization of the baro- 
metric law. 

If the plasma is not collision dominated, 
randomization of electron velocities may 
still occur through the enhanced levels of 
fluctuating fields in the plasma, such as 
occur for electron two-stream instabili- 
ties, or electron-ion instabilities of the 
lon-acoustic or Bunemann ti'pe.(9,9) These 
mechanisms are probably effective in coup- 
ling neutralizer electrons into the bulk 
plasma and in ecualizing electron and ion 
mean drift velocities. They are_often 
approximated by introducing an effective 
collision frequency, j, in place or 

The determination of electron tempera- 
tures in tne plasma requires considera- 
tion of the energy balance equation, equa- 
tion (5) . .Making the sane approximations 
in the equation expressing conservation of 
energy that were made in tne momentum 
equation, yields 
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7-F 


(13) 


» qn£-V + a -V 
a a 


+ Q 


ai 


with 

F - I ?V q (14) 

Here 3/2 pV is the enthalpy flux of the 
drifting electrons, q the nacroscopic 
heat flux, and R*V is related to the af- 
fective joule heating associated with the 
relative notion of electrons and ions. 

The quantity R appears also in the elec- 
tron momemtum equation: for a plasma con- 
trolled by collective affects it should 
be approximated in the energy equation 
in the same manner as in the momentum 
equation. The heat flux, q, contains new 
features. Classically, 5 contains two 
terms; one proportional to the relative 
drift velocity between electrons and ions, 
and the other proportional to the gradient 
of electron temperature. 

For the initial calculations, we ignore 
the drift contributions to the energy flux, 
the electron- ion heating Qei» and assume 
that the heat flux is proportional to the 
temperature gradient. The energy balance 
equation thus assumes the simple form 

7*<70 + mnv(Vg-V^)“ => 0 (IS) 
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4. Ion Snqine Neutralizer Code 

The basic physics of an ion engine 
neutralizer model was presented in the 
preceding sections. This physical model 
has been incorporated into a two-dimension- 
al (R-Z) computer coda, which is described 
below. A seunple calculation of neutrali- 
zation in a t.hruster similar to SERT II 
IS discussed below. Results are given 
for space dependent electric potentials, 
electron te.mperatures , and current densi- 
ties. 

Coda Description 

The icn t.hruster model has been incor- 
porated into a two-dimensional (R-Z) com- 
puter coda following the bloc)c diagraua 
shown i.n Figure 1. For t.his initial 
version t.he ion currents anc densities 
ware assum.ed known . In a later version. 

It would be possible to allow a multi- 
component ion composition to be determined 
self-consistently with the temperature and 
potential. The code operates entirely in 
.MKS units. 


Fig. 1 Bloc.k diagram for ion engine 
neutralization code. 

At present, the code assumes ion veloc- 
ities everywhere and ion currents at the 
lapum boundary to oe )cnown. The code then 
calculates plasma densities such that 
7-(ny) = 0 13 numerically satisfied. Typi- 
cally, ion velocities are taken to be either 
purely axial or to be radial from a point 
source on Che axis exterior to the mesh. 

As the code requires ton-zero plasma den- 
sity everywhere, a background density of 
"slow" ions may be added. It should be 
possible to handle multiple ion species 
with interconversion fairly easily. 

The neutralizer is assumed to be a ring 
at specified distance from the axis, emit- 
ting a currant of electrons equal to the 
ion beam current. The net current i.n the 


plasma 

IS given by 

(16) 

2 = 

nq(Vi - v^) 

= 

a(-7 3 - 2 Tp) , 

(17) 


The code has been run interactively, 
with all relevant information on disk file 
As long as previous information exists on 
disk, the program may be entered from t.na 
two noted entry points as well as t.he be- 
gi.nning. For developmental purposes, it 
was found convenient to "hard-wire" many 
features of a particular problem into the 
code, while others are prompted for .nput. 
.A flexible on-lj.ne graphics program, whic.h 
plots information on the disk file, .nas 
also been developed. 


where now p = nkT is the electron pressure. 
(For r— ” and qkT— ’ (constant) we find 
j = -In n.l The code determines electro- 
static potentials by solving "'•j = 0. (See 
Appendix A.) It 13 necessary to iterate 
between this equacion and the temperature 
equation (equation 15) , since the pressure 
IS a function of temperature. 

The plasma temperature satisfies t.he 
equation 
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7- (-<79) 



7 


(13) 


where < is the theraial conductivity and 
the right-hand-side represents the ohmic 
generation of heat. For this preliminary 
version, convective heat transport has 
been neglected. On the various boundary 
regions, either isot.hermal or insulating 
boimdary conditions nay be specified. 

Since, in practice, we take < to have a 
power law dependence on 9 , < = nic'9'^~^, 
tne equation actually solved is 

2 

-^•<'7(9''*) » ^ . (19) 


For convenience, the transport coeffi- 
cients a and c' are calculated by a single 
isolated subroutine. The conductivity a 
may depend on both density and temperature, 
and <' on density only. The present ver- 
sion assumes a relaxation race proportional 
to the plasma frequency: 


1/2 e‘ 1 
" m 3.98a 


( 20 ) 


where the parameter a is ta)cen to be 
O.Sl. 3y the classical Weideman-Franz 
law. 


the beam from t.he side. The ion currant is 
SO percent neutralized at ^13 cm downstream 
from the thruster. 



Fig. 2 Plasma density for SERT II ion 
thruster model. 



If we measure temperature in eV, )< « q, 
so that <' = 3/4 c. 

5. Computational Results 

A calculation was performed for neu- 
tralization of a 0.23 ampere, purely axial 
beam of 3 ksV mercury ions in a constant 
density background plasma. The beam had 
a radius of ^7 cm, and a ring neutralizer 
was placed at a 17 cm radius. The given 
plasma density (Figure 2) had a peak of 
16.1 X 101“* m"^, and an ambient density 
of 1.0 X 10-^ ■* . These conditions ap- 

proximate those occurri.ng in the plasma 
produced by the SERT II t.hruster. The 
temperature profile (Figure 3) was calcu- 
lated witn insulating boundary conditions 
at the thruster. The maximum temperature 
occurs at the beam entrance , where the 
heat generation is greatest. The peak 
temperature was S eV, compared with a 1 eV 
assumed background. The electrostatic 
potentials are shown in Figures 4 and 5. 
Strong potential variations across the 
beam edge and a potential dip near the 
neutralizer are calculated features i.n 
qualitative agreement with experimental 
results (Figure 5) (ID While the strong 
edge fields conform approximately to the 
barometric law at t.he local temperature, 
they deviate substantially from the results 
that would be obtained wit.n v << jp or by 
using the baror.etric law as the point of 
departure (t.hat is, a zero resistivity, 
isothermal calculation) . Current-vector 
plots (Figures 6 and 7) indicate that t.he 
beam is neutralized by electrons entering 
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Fig. 3 Temperature profile for SERT II 
ion thruster model. 
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P0TENT1RL3 (eV) 


LOO NET J 


R-JVXIS (m) 

Pill MAX - 21.40 

rig. 4 Electrostatic potentials for 
SEST II ion thruster 'nodal. 


R-AXIS (m) 

Current vectors for SERT II _on 
thruster •nodel. The length of each 
arrow is wea/cly dependent on cur- 
rent. 

NET R H J 


POTENTIALS :• a.oo 



R-«IA (n) 

Electrostatic potentials for 
SERT II ion thruster model. 


■i 11/ 

•I// — 

3 

ti M U Zl Zt 

R-AXIS (m) 

Fig. 7 Current vectors for SERT II ion 

thruster model. The length of each 
arrow is proportional to radius 
tunes current density. 

6. Sunmarv anc Conclusions 

The consequences of nhe assumption that 
the electron gas .near an ion thrister be- 
haves as a resistive fluid have been ex- 
aminea. Theoretical results obtained here 
are in qualitative agreement with experi- 
mental observations. Such agreement indi- 
cates that the orooerties of thruster 





generated plasmas can be described by 
fluid equations having a classical form 
but with an effective collision frequency 
near much in excess of the classical 
value iOr pair (electron- ion) collisions. 

Further worh should be performed to 
test qualitative and queuititative predic- 
tive capabilities of fluid models of 
thruster plasmas, and to better understand 
the relationship between ersatz fluid 
models of the type invoked here and the 
underlying plasma physical mechanisms 
embodied in the collisionless Vlasov 
equation. 

Appendix A. Mathematical Considerations 
' ' on the Variational Formulation of' 
Poisson's Equation 


where 

-/ d^r J(r) [7N,(r)l • (TN, (r)l (A. 5) 

Sj =•/ d^r Nj(r) S(r) , (A. 6) 

1 

Sj_(t) are nodal interpolation functions 
and'pj_ are node values for the unknown 
field. 

Since W IS symmetric, these aquations 
are = 

LWij 4>j + Si = 0 (A. 7) 

Using (A. 5) we find 


In our theory of ion engine neutraliza- 
tion, It IS .necessary to solve equations 
of the form 

7- [o(r)7 z] = S(r) (A.l) 

subject to fixed-value ooundary conditions 
at some nodes ^u^d normal-gradient boundary 
conditions at others. Me need to show 
that aquation (.A.l) is exactly the equa- 
tion equivalent to the variational formula- 
tion, and that the normal-gradient bound- 
ary conditions are equivalent to a surface 
c.narge. 

Theorem 1 

Minimization of 

T I 7(r) , ) 

j^d-^r ’ 7.*!“^ t t(r)S(r) j (A. 2) 

IS equivalent to equation (A.l) . 


2 */ d^r(7N (r) I • [T(r)": (r) 1 

J ‘ (A. 3) 

where 

T(r) * 2 M (r) 

1 * ' “ 

The integrand of (A. 3) can be rewrit- 
ten as 

7* tc(r)Nj_(r)7j(r) ] 

- Mj_(r)7* [3(r)7l(r) ] (A, 9) 

Writing B(q) » c(r)'3~(r), and using the 
divergence theorem, aquation ^A.7) be- 
comes 

/d^r M^(r) D(r)*n -J^d^r M^_(r) 7.o(r) 


Proof 

Minimization of fd?zHz ,~i ,s) requires 
31(t,"j,r) 3(.(:,7c,r) 

' ^ =• 0 (A- 

Taxing L as t.he i.ntegrand of (A. 2), equa- 
tion (A. 3) yields exactly equation (A.l). 

Theorem 2 


-t- » 0 (A. 10) 

where the first i.ntegral runs over the 
surface of volume 7. By analogy with 
electrostatics, the second and tnird 
terms are the "volume charge" associated 
with noae i, wnila the first, involving 
the surface-normal-gradient of the poten- 
tial, IS non-zero only on surface nodes 
and may be compensated by a corresponding 
surface-charge density. 


In the finite-element formulation for 
the minimization of (A. 2) , normal-gradient 
boundary conditions are equivalent to sur- 
face charge. 

Proof 


In the finite-element formulation, 
.minimization of (A. 2) leads to the equa- 
tions 



z 

3K 


'-"jk ^k " 




0 


(A. 4) 


Corollary 


Consider the equation 

7-J =• 0 (A. 11a) 

where 

J = c(7i-g(r)7f) , (A. 11b) 

and f(r), g(r) are known functions. If 
the inhomogeneous term is evaluated using 
equation (A. 3), equation (A. 11), wit.h no 
further qualifications, will give zero- 
nomal-current bouncary conditions when 
solved by the finite-element method. 
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rurthernore, non-zero ^noraal current 
boundary conditions, n J 3 (r), can be in- 
voked by adding an inhomogeneous term to 
the surface node aquations. 

Proof 

Substituting (A. 11) into (A. 10), the 
finite element equations are 

Td^r N (r) J(r)*n - N (r) » 0 

S'l '■ --- 1 - 

(A.12) 

As the second term is to vanish (within 
the finite-element approximation) , it fol- 
lows that the first terra is equal to the 
right-hand-side. Thus, replacing the 
right-hand-side with 

f d^r M^(r) Jg(r) (A. 13) 

S 

will produce the specified boundary condi- 
tions. 
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Systems , Scxence amd Software 
La Jolla, Calxfomia 92038 


Abstract 


Solar electric propulsion is a leading 
candidate for many upcoming space missions. 
Under many circumstances plasma produced by 
charge-exchange reactions within the ion 
beam dominates the ambient environment near 
the spacecraft. The calculations presented 
here contain a predictive hydrodynamic 
model for the charge-exchange plasma ex- 
pansion, and a fully three-dimensional 
model for the structure of the plasma 
sheath around the solar array wing. Re- 
sults of calculations for several configu- 
rations and voltage levels indicate that 
with kilovolt biases power losses of 'vlO 
percent or more are likely, even with only 
one engine in operation , and that amelio- 
rative measures should focus on the inboard 
portion of the solar arrays . 


E 



„(P) 

a 


Ar 

Az 

9 

P 

A 


Nomenclature 


= electric field 
= radial mass flux 

= axial mass flux 

= acceleration on 

= Boltzman's constant 
= ion mass 
= electron density 
* ambient plasma density 

= magnitude of electron charge 
= radial coordinate 


radial location of velocity vector 


V 


13 


= radial location of ion density p 

13 


k'"^ surface boundary of n 


13 


= electron temperature 
= ion velocity 
= axial coordinate 


= axial location of velocity vector 

^13 

= radial location of ion density p 


13 

angle from the thruster beaun direc- 
tion 

radial mesh spacing 
axial mesh spacing 
kT 


= ion density 

= electric potential 

= volume of cell associated with V 

13 


1 . Introduction 


Solar electric propulsion is a leading 
candidate for lifting large structures 
from shuttle orbit to geostationary alti- 
tude. Previous studies demonstrated 


♦Program Manager, **Senior Research Scien- 
tist, tResearch Scientist 


that plasma produced by charge-exchange 
reactions within the ion beam may dominate 
the ambient environment near the space- 
craft. Currents flowing through this 
plasma between the engine neutralizers and 
the solar arrays are a drain on the power 
systems. If the losses become too large, 
they may have a substantial mission impact. 

Simple calculations have been performed 
to estimate parasitic currents flowing 
through the charge-exchange plasma. 

While the potential seriousness of the 
interaction was identified, the' ad hoc 
nature of the previous work made clear the 
need for a more accurate treatment of the 
expansion of the charge-exchange plasma 
and the resultant solar array power losses. 
The calculations presented here are an im- 
provement over previous work in that they 
contain a predictive model for the oharge- 
exchainge plasma expansion, and a fully 
three-dimensional model for the structure 
of the plasma sheath around the solar array 
wing. 

The generation mechanism of the charge- 
exchange plasma is clearly understood. 3 
However, until recently there were no 
predictive models of the expansion dynamics 
of these slow moving 10ns. This work, as 
well as that of Robinson, Kaufman and 
Winder, 4 describe the 10ns as streaming 
under the electric field created by baro- 
metric law behavior of the plasma elec- 
trons . The model described here treats 
the 10ns numerically as a cold fluid while 
Reference 4 calculates representative ion 
trajectories by particle pushing. Both 
calculations are presently done assuming 
axial symmetry. 


The collection of 10ns and electrons by 
solar arrays has been the subject of much 
recent interest. Experiments have been 
performed on large array- like objects in 
the vacuum tank at NAS A/ JSC. 5 a three- 
dimensional computer code, which simulates 
plasma interactions in Low Earth Orbit 
(LEO) was developed by the authors and 
provided qualitative and quantatitive 
agreement with experiment. ® This pro- 
grcun, NASCAR /LEO, in a somewhat modified 
form IS used in this study to describe 
the charge-exchange plasma-solar array 
interactions. 


In this paper we present a new model 
for the expansion of the charge-exchange 
plasma from an ion thruster. We then use 
this predicted charge-exchange plasma as 
the environment surrounding solar arrays 
in the LEO plasma interaction model. From 
this we obtain improved estimates of power 
losses for a variety of solar array con- 
figurations . 


Copyriglit © Amencan Institoie of Aeronautics and 
Astronautics. Inc., 19»1 Ail rights reserred. 
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2 ■ Theory 


3 . Numerxcal Methods 


For this study we have considered a 
single 30 cm dicuneter mercury thruster 
producing 2 A of beam current and 25 mA of 
charge-exchange current- ^ The charge- 
exchange ions are emitted radially, with 
energies of 5-10 eV, within a downstream 
distance equal to one diameter of the beam. 

14 -3 

This gives a density of ”',2.5 x 10 m 
at the beam edge, which exceeds the Low 
Earth Orbit (LEO) plasma density of 

10^°-10^^ m“^. Thus the charge-exchange 
plasma will dominate the ambient to a 
distance of several meters from the 
thruster in LEO, and over the entire space- 
craft for substantially more tenuous en- 
vironments . 

The conditions in this plasma are long 
collision lengths {''-lO^ m) and short Debye 

lengths (''-10~^ m) . In this regime we 
assume that the electron gas adjusts it- 
self to maintain isothermal, quasi-neutral 
conditions. This implies a barometric law 
potential , 

(J> = Qln{n/n^) (1) 

where 9 is the electron temperature (eV) , 
n the local plasma density, amd n^^ the am- 
bient plasma density. The barometric law 
IS not an essential feature of the model. 

In the future we expect to remove the iso- 
thermal restriction and utilize fluid-like 
equations to describe the electrons. ^ 

The barometric law potential causes expan- 
sion of the ion cloud as it emerges from 
the beam edge. Far downstream the electro- 
Static force become small r so that the ion 

2 

density takes the form f(a)/r , where a is 
the polar angle relative to the beam and r 
the distance from the engine. 

The ion gas is modeled hydrodynamically , 
i.e., we assume the ion density and velocity 
to be a well-defined function of position, 
and the ion thermal motion to be unimport- 
ant. This representation is chosen for 
ease of generalization, as opposed to using 
particle pushing techniques which require 
extensive computer time for three-dimen- 
sional applications. The motion then 
satisfies the equations of continuity of 
mass and momentim: 


In order to suppress spurious numerical 
oscillations we solved equations (2) amd 
(3) using upwind differencing on a stag- 
gered mesh. ® The computational mesh is 
shown in Figure 1. The thruster beam 
starts on the z =• 0 plame and is pointed 
along the positive z aocis. Charge-exchange 
ions are produced downstream within the 
thruster beam, and the charge-exchange ions 
are expelled from the main beam radially by 
electric fields within the beam. The phy- 
sical space IS located within the solid 
boundary? mesh points outside the physical 
space are computationally convenient and 
serve to maintain boundary conditions. 
Velocities are defined at the points indi- 
cated by crosses, which form an evenly 
spaced r-z mesh: 



Fig. 1 Computational space, showing points 
for definition of velocity (x), 
density (•) , and mass flux (-»•} . 



-7- (pv) 


( 2 ) 



r^ + 1 Ar 


(4a) 


9t 


(pv^) 


-y.(v^Pv) ^ 


(3) 



z + 1 Az 
o 


(4b) 


where subscript a denotes a Cartesian com- 
ponent, and E IS the electric field. The 
challenge is to develop numerical methods 
capable of finding the steady-state solu- 
tion to equations (1-3) in the R-Z 
(cylindrical) geometry appropriate to the 
problem. 


Densities are defined at the centroids of 
the quadrilaterals formed by the crosses: 


,(P)2 

11 

(,,V.2 , ,,vu) 

(5a) 

(P) _ 
3 


>r’ * 

(5b) 
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Mass fluxes, f, are defined at the arrows 
in a fully "upwind" sense: 




■ K'V'r’) 


k = 


1-1 if Vj. ^ 0 

1 if < 0 




j 3-1 If ^2 i 0 

(3 


If < 0 


With the above definitions, it is straight- 
forward to timestep the integral form of 
equation (2) , using as control volumes the 
rectangles formed by the four crosses 
(velocity definition points) surrounding 
a density definition point. 

Equation (3) is then used to construct 
dv/dt. We use a finite difference method 
to evaluate the divergence term as follows: 

a. Define mass fluxes flowing between 
the crosses in Figure 1, i.e., normal to 
the arrows. This is done by taking 
weighted averages of nearby already de- 
fined fluxes: 




1+1 3 

E E 

k=l S,=3-l 




(7) 


b. Define densities associated with 
velocity definition points: 




(6a) 

i 

E 



jC=1-1 

i.j-i 


These densities are needed at both their 

(v)\ 

previous and current (just updated) values. 

3 / 

G. Define 

velocities at the arrows: 

(6b) 

(^1+1 

,<V,) 


S 



+ 


(6c) 

r3+i 

«r) 


\ 

- »r’) 

"3 ) 

S 

(6d) 

+ 



d. Determine the change in velocity by 


^13 3t 


.(v) ,(v)\ 


- - * ”4 


where is the volume 


< 4 < 

3 aj^e the surfaces bounding that volume, 

/ (v) (v)\ 

and is the force at Ir^ ,z^ j . 

Specifically, using 

we obtain the following expression for the 
divergence of the momentxim flux integrated 
over the volume: 


F 

3 -13 


(P) 




3+1 


E E4“’«.(4“>-r’) 

k=i-l 1=3 


(8) 


- ^ =k "k-!k :fk 
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The velocities are advanced using 


Vi^ (t) = Pj_j(t-dt) v^j (t-dt)/p^^ (t) 

+ (q/ra) E dt 

■ n P^^(t) 2 ^kf^Tc'fk ^k) • 

ij 13 ^ \ / 


These formulae are accurate to first order 
in dt and second order in {^r,Az) . 

e. Additionally, provision is made for 
optional velocity smoothing at each time- 
step. 

The problem is run until there is no 
variation of velocities in time. This 
empirical steady state is generally 
achieved in a few times the time it takes 
for a fluid particle to traverse the mesh. 

4. Results 

We have applied the above model to cal- 
culate the expansion of an ion charge- 
exchange plume with initial conditions 
similar to a case measured by Hughes Re- 
search Laboratory. 3 The initial condi- 
tions were: 

Emission energy (radial) = 10 eV. 

Initial radius = .15 m. 

Current density q 0 zj , z j 

( 0 z<0 

3 ( 

(.123 exp(-z/.22) z>0 

(total current » 25 mA) 

Electron temperature 9 * 1 eV. 

Mass (Hg ion) = 3.34 x lo"^^ kg. 

12 -3 

Ambient density n^ = 10 m . 

It follows that 



= 3095 

m/sec 


= 2.48 

X 10^'* 


*max = = 5-5 V. 


The calculation was done in three 
phases: 


a. 

.15 

m 

< 

r 

< 

1,15 m; Ar = 0.05 m 


-1.0 

m 

< 

2 

< 

1,0 m; Az = 0.10 m 

b. 

1.0 

m 

< 

r 

< 

5.0 m; Ar = 0.2 m 


-2.0 

m 

< 

z 

< 

2.0 m; Az = 0.2 m 

c« 

2.0 

m 

< 

r 

< 

12 m; Ar = 0.5 m 


-4.0 

m 

< 

2 

< 

4m; Az » 0.4 m 


Initial conditions for the second and 
third phases were obtained by interpolat- 
ing data from the previous phase and re- 
normalizing to retain the correct total 
current. At each phase, the calculation 
was carried out until steady-state den- 
sities and velocities were reached. 

The results are shown in Figures 2 and 
3. It IS seen that the initially asym- 
metric expansion (Figure 2) becomes 
roughly symmetric by a radius of 'vl m from 
the beeim. This is seen in Figure 3a, in 
2 2 

which spherical (r + z ) scaling maintains 
constant arrow length, and in Figure 3b, 
which indicates radial contours of equal 
2 2 

(r + z ) times the density. The density 

(m"^) beyond 1 m from the engine is 
reasonably approximated by 

p = 10^^ (sina)^^/(r^ + z^) 

where a is the angle from the beam direc- 
tion. A closer inspection of Figure 3b 
indicates that the plume extends further 
upstream than downstream. This is attri- 
butable to "pressure blowoff" at the up- 
stream plume edge near the engine where 
the density gradient is high. 

The plasma parameters determined from 
the above model are used as input to a 
fully three-dimensional computer program 
designed to predict current collection by 
high voltages in low temperature, short 
Debye length plasmas. ^<3 This model uses 
an analytic , nonlinear space charge formu- 
lation, correct in both Debye screening 
and thin sheath limits , to determine the 
electrostatic potential and the boundary 
of the plasma sheath. The model allows 
the plasma temperature and density to vary 
in space. Figure 4 gives a sample of 
electrostatic potential contours near a 
solar array wing, illustrating the asym- 
metry caused by the charge-exchange plasma 
being predominantly on one side of the 
wing. By tracking electrons inward from 
the sheath boundary, good parasitic cur- 
rent estimates are obtained. Iterating on 
the two stages of the calculation allows 
non-local effects to be included. Addi- 
tionally, a high resolution capability is 
available to compute the current distri- 
bution over a complex pattern of solar 
cells. 

Table 1 gives sample results for an 
8 m X 30 m, 25 kw solar array for several 
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. 2 Cylindrically scaled current den- 
sity plots for ion plume flow (a) 
near the beam, and (b) in the 


"spherical expansion" region. 





R AXIS (HETXRS) 


Fig. 3 Spherically scaled plots of (a) 
current density auid (b) particle 
density contours. 



Fig. 4 Potential contours about a 2 kV 

solar array in the presence of the 
charge-exchange plume. Note that 
the potential is more shielded 
above the wing (where the charge- 
exchange plasma is located) them 
below. The contour differences 


are 100 volts. 


orientations and configurations. In these 
calculations, the spacecraft body was held 
at plasma ground and the solar aurray was 
divided into three equal sections, each of 
which was positively biased. Parasitic 
currents to both sides of each 10 m sec- 
tion of the wing were separately monitored. 
Several trends are apparent; (1) The in- 
board section of the wing draws most of 
the current, even though it is at a low 

voltage. This is because of the r~^ ex- 
pansion of the charge-exchange plasma. 

(2) The current to the outbocurd section 
IS similar to the center section current, 
due to the large end effect in the tenuous 
plasma. (3) When the beam is in the plane 
of the panel, an increased loss is caused 
by the array's intersecting the charge- 
exchange pancake. It is apparent from 
these results that with kilovolt biases 
power losses of '^10 percent or more are 
likely, even with only one engine in opera- 
tion, and that ameliorative measures should 
focus on the inboard portion of the solar 
arrays . 

5. Discussion 

As large, high powered solar electric 
propulsion vehicles come closer to reality. 
It becomes more important to place realistic 
bounds on the plasma array interactions. 
Previous studies ^ad assumed simple 
models of the charge-exchange expansion. 

The expansion model is so important in 
estimating power losses that the authors ' 
previous work 2 parameterized the results 
in terms of an undetermined expansion 
angle. The calculations presented here 
are a vast improvement, inasmuch as the 
expansion of the charge exchange plasma is 
determined on a reasonable theoretical 
basis and the plasma collection is fully 
three-dimensional. 


However, certain other aspects of the 
model presented here are just as primitive 
as in the earlier studies . In particular 
the isothermal barometric law description 
of the electrons, the lack of self- 
consistency between expansion and sheath 
models, and the absence of any turbulent 
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Table 1 Power losses for 3 m x 30 m solar array for 
one thruster and various configurations 


Configuration 

Inboard 

Voltage 

Center 

Outboard 

Current ( one 
Inboard Center 

side) 

Outboard 

Power Loss 
(per wing) 

a,d 

1000 

2000 

3000 

1.4 

0.6 

0.6 

4.4 

b,d 

2000 

2000 

2000 

4.7 

0.6 

0.5 

12.0 

b,d 

1000 

2000 

3000 

3.7 

0.6 

0.5 

6.4 

a,c 

1000 

2000 

3000 

1.1 

0.3 

0.4 

5.9 

a.,c 

100 

1000 

2000 

0.3 

0.3 

0.4 

2.3 

a,d 

200 

200 

200 

1.4 

0.3 

0.2 

0.38 

a,d 

400 

400 

400 

1.5 

0.3 

0.3 

0.85 


a - Beam in plane of panel 

b - Beam normal to plane of panel All voltages in volts, currents in amps, 

c - Both sides biased euid power losses in kilowatts, 

d - Back of panel grounded 


heating mechanism during electron col- 
lection place severe restrictions on the 
accuracy of these results. For mission 
analysis, there also remains unanswered 
the question of charge-exchange plasma 
generation and expansion in a multiple 
thruster configuration. These are some 
of the technical questions which must be 
addressed if we are to provide an accurate 
assessment of parasitic current losses due 
to ion propulsion generated plasmas. 

References 

1. Kaufman, H. R. , "Interaction of a 
Solar Array with an Ion Thruster Due 
to the Charge-Exchange Plasma," 

NASA CR-135099, 1976. 

2. Parks, D. E. and I. Katz, "Spacecraft- 
Generated Plasma Interactions with 
High-Voltage Solar Array," J. Space- 
craft and Rockets, 16, 259, 1979. 


7. Parks, D. E., M. J. Mandell and 

I. Katz, "Fluid Model of Neutralized 
Ion Beams," AIAA-81-0141, AIAA 19th 
Aerospace Sciences Meeting, St. Louis, 
MO, January 12-15, 1981. 

3. Richtmyer, R. D. and K. W. Morton, 
Difference Methods for Initial-Value 
Problems , Interscience Publishers , 

New fork. Second Edition. 

9. Msuidell, M. J. , I. Katz, P. G. Steen 
and G. W. Schnuelle, "The Effect of 
Solar Array Voltage Patterns on Plasma 
Power Losses," IEEE Trans. Nuc. Sci. , 
Vol. NS-27, No. 6, p. 1797, December 
1980. 

Aclcnowledgments 

This work supported by NASA/Lewis Re- 
search Center, Cleveland, OH under Con- 
tract NAS3-21762. 


3. Poeschel, R. L. , E. I. Hawthorne, 
et al. , "Extended Performance Solar 
Electric Propulsion Thrust System 
Study," NASA CR135281, 1977. 


4. Robinson, R. S., H. R. Kaufman and 
D. R. Winder, "Simulation of Charge- 
Exchange Plasma Propagation Near an 
Ion Thruster Propelled Spacecraft," 
AIAA-81-0744, AIAA/JSASS/DGLR 15th 
International Electric Propulsion 
Conference, Las Vegas, NV, April 21-23, 
1981. 


5. McCoy, J. E. and A. Konradi, "Sheath 
Effects Observed on a 10 meter High 
Voltage Panel in Simulated Low Earth 
Orbit Plasma," Spacecraft Charging 
Technology-1978, 315, NASA Conference 
Publication 2071, AFGL-TR-79-0082 , 
1979. 


6. Katz, I., M. J. Mandell, G. W. 

Schnuelle, D. E. Parks and P. G. Steen, 
"Plasma Collection by High Voltage 
Spacecraft at Low Earth Orbit," J. of 
Spacecraft and Rockets, Vol. 18, No. 1, 
p. 79, January-February 1981. 


84 


6 


APPENDIX C 


NEUTRAL DENSITY FREE EXPANSION MODEL 

Problem ; A Maxwellian spectrum of velocities is emitted 
uniformly from a disk of radius a located on the plane 
z = 0. Find the density of particles throughout space. 

Geometry: 



Analysis : On a trajectory 

f (p,z,V|^,v^) = f^ (p^,o,V|^,v^) 

where f is the distribution function at r. 
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since 


-3/2 


= 2n^(TTV^) 


exp 


2 2 
Vf + V 

I z 


V! 


p < a, V > 0 


= 0 


p > a or V <0 
2 


n = density 
o 


2 = 0, p < a 


= 2 kT/m 


: and are connected by a trajectory 


p - Pq = 2/v^ 


V 


d^v, dv^ = -f d^ 


-3/2 


n = 


2n^(7rVj) 




(P-Pq) 


+ 1 


dv 


= 2n^(TTVj) 


-3/2 


= 2n^(7TV;) 


-3/2 




exp{- 


^ \ 2 2 
(p-p^) + 2 


2 2 , 

X > X dx 


/tt 

4 


(P-Pq)^ + 2^ 


.-3/2 


a 27 t 


•'o •'o 


PodPod4> 


2 2 2 

p + pQ + 2 - 2pp^ COS(j) 


-3/2 


a 2TT 

W (- h)f f ‘’o'SOo'^* 

*'o •'o 


2 2 2 

P + P +2 - 2pp cos4> 

o o 


- 1/2 
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Observations : 


1. For p < a, the integral is singular for z ^ 0; 
n(p,0) = n^ for p < 0, n(p,0) = 0 p > a. 

2. The density field is formally equivalent to the E 

z 

field of a disk with a uniform distribution of surface 
charges. 

2 2 2 2 

3. Asymptotically, (p +z =n >>a) 

n z/r^ 


Integration : The integral over p^ can be performed 

analytically 


2ir 

" = 2 ? 

•'A 


-Ir 


di^ (z^+p^sin^({) ) 


r - 


r -ap coscj) 


2 2 

(r +a -2apcos4>) 


1/2 


where 


/ 2 2 , 

r = (p + z ) 


1/2 
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